{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Confidence interval of survival data\n", "\n", "## Introduction\n", "\n", "**Survival analysis** is a branch of statistics that deals with the analysis of time-to-event data, often referred to as 'survival data'.\n", "\n", "The **event** can be anything from the death of an organism (hence **survival**), to the failure of a mechanical system, to the recovery of a sick patient.\n", "\n", "A **confidence interval** in survival analysis gives a range of values, derived from the data, that is likely to contain the true value of an unknown population parameter. In the context of survival analysis, a confidence interval might be used to describe the uncertainty around a survival time estimate.\n", "\n", "## Survival data\n", "\n", "### The `lifelines` library\n", "\n", "In Python, the [`lifelines` library](https://lifelines.readthedocs.io/en/stable/index.html) is commonly used for survival analysis. It provides a simple and intuitive API for fitting survival models, and also includes functions for visualizing survival data. Here's a basic example of how we might calculate a confidence interval for a survival function estimate using the `KaplanMeierFitter` class in lifelines:" ] }, { "cell_type": "code", "execution_count": 42, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "lifelines version: 0.29.0\n" ] } ], "source": [ "import pandas as pd\n", "import matplotlib.pyplot as plt\n", "import lifelines\n", "\n", "print('lifelines version:', lifelines.__version__)" ] }, { "cell_type": "code", "execution_count": 43, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
YEARSEVENT
16.540
36.170
63.670
04.071
21.391
45.891
54.761
\n", "
" ], "text/plain": [ " YEARS EVENT\n", "1 6.54 0\n", "3 6.17 0\n", "6 3.67 0\n", "0 4.07 1\n", "2 1.39 1\n", "4 5.89 1\n", "5 4.76 1" ] }, "execution_count": 43, "metadata": {}, "output_type": "execute_result" } ], "source": [ "from lifelines import KaplanMeierFitter\n", "\n", "# Assume we have some survival data in the following two lists:\n", "# `durations` is a list of durations until the event or censorship\n", "# `event_observed` is a binary list where 1 indicates the event, \n", "# i.e., death, occurred, and 0 indicates censorship\n", "\n", "# example of the survival data from pages 47-48 of the book Intuitive Biostatistics\n", "durations = [4.07, 6.54, 1.39, 6.17, 5.89, 4.76, 3.67]\n", "event_observed = [1 , 0 , 1 , 0 , 1 , 1 , 0 ]\n", "\n", "# we create a dataframe with those data\n", "data = pd.DataFrame(\n", " {\n", " 'YEARS': durations,\n", " 'EVENT': event_observed,\n", " }\n", ")\n", "\n", "# sorted durations\n", "data.sort_values(by='EVENT')" ] }, { "cell_type": "code", "execution_count": 44, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 44, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Create a fitter instance\n", "kmf = KaplanMeierFitter()\n", "\n", "# Fit the data to the model\n", "kmf.fit(\n", " durations=durations,\n", " event_observed=event_observed,\n", ")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The `kmf._median` attribute gives us the **median survival time**. In the context of survival analysis, the median survival time is the smallest time at which the survival probability drops to 0.5 or below. In other words, it's the time by which half of the population has experienced the event of interest." ] }, { "cell_type": "code", "execution_count": 45, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Median survival time: 5.89\n" ] } ], "source": [ "print(\"Median survival time:\", kmf.median_survival_time_)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Please note that the median survival time might be `inf` (infinity) for some datasets. This happens when the survival function never goes below 0.5, which means that more than half of the population survives till the end of the observation period. In such cases, the median survival time is undefined and is conventionally considered to be infinite." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The survival function estimate and its confidence intervals can be accessed with:" ] }, { "cell_type": "code", "execution_count": 46, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
KM_estimate
timeline
0.001.000000
1.390.857143
3.670.857143
4.070.685714
4.760.514286
5.890.342857
6.170.342857
6.540.342857
\n", "
" ], "text/plain": [ " KM_estimate\n", "timeline \n", "0.00 1.000000\n", "1.39 0.857143\n", "3.67 0.857143\n", "4.07 0.685714\n", "4.76 0.514286\n", "5.89 0.342857\n", "6.17 0.342857\n", "6.54 0.342857" ] }, "execution_count": 46, "metadata": {}, "output_type": "execute_result" } ], "source": [ "kmf.survival_function_" ] }, { "cell_type": "code", "execution_count": 47, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
KM_estimate_lower_0.95KM_estimate_upper_0.95
0.001.0000001.000000
1.390.3340540.978561
3.670.3340540.978561
4.070.2127970.912112
4.760.1177600.813249
5.890.0481080.685484
6.170.0481080.685484
6.540.0481080.685484
\n", "
" ], "text/plain": [ " KM_estimate_lower_0.95 KM_estimate_upper_0.95\n", "0.00 1.000000 1.000000\n", "1.39 0.334054 0.978561\n", "3.67 0.334054 0.978561\n", "4.07 0.212797 0.912112\n", "4.76 0.117760 0.813249\n", "5.89 0.048108 0.685484\n", "6.17 0.048108 0.685484\n", "6.54 0.048108 0.685484" ] }, "execution_count": 47, "metadata": {}, "output_type": "execute_result" } ], "source": [ "kmf.confidence_interval_" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This gives us the survival function estimates along with their confidence intervals. The first table corresponds to the \"% survival\", and the second table to the \"lower limit\" and \"upper limit\" in the table 5.3 of the book.\n", "\n", "Remember, the interpretation of confidence intervals in the context of survival data is the same as in other areas of statistics. A 95% confidence interval means that if we were to repeat our study many times, we would expect the true survival function to fall within our estimated interval in 95% of studies." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The `lifelines` API provides many other utilities, for example the function called `survival_table_from_events` is used to create a survival table given some durations and censoring vectors. The survival table provides a summary of the number of individuals at risk, the number of events, and the number of censored observations at each unique time point in the data." ] }, { "cell_type": "code", "execution_count": 48, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " removed observed censored entrance at_risk\n", "event_at \n", "0.00 0 0 0 7 7\n", "1.39 1 1 0 0 7\n", "3.67 1 0 1 0 6\n", "4.07 1 1 0 0 5\n", "4.76 1 1 0 0 4\n", "5.89 1 1 0 0 3\n", "6.17 1 0 1 0 2\n", "6.54 1 0 1 0 1\n" ] } ], "source": [ "from lifelines.utils import survival_table_from_events\n", "\n", "# Create the survival table\n", "table = survival_table_from_events(durations, event_observed)\n", "\n", "print(table)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This prints a DataFrame that includes the following columns:\n", "\n", "- _removed_ is the number of individuals that had an event or were censored at each time point.\n", "- _observed_ is the number of individuals that had an event at each time point.\n", "- _censored_ is the number of individuals that were censored at each time point.\n", "- _entrance_ is the number of individuals that entered the risk set at each time point. This is usually 0 for all time points except the first one.\n", "- _at_risk_ is the number of individuals that are at risk of having an event at each time point.\n", "\n", "And `lifelines.utils.survival_events_from_table` is the inverse of the function survival_table_from_events:" ] }, { "cell_type": "code", "execution_count": 49, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Durations of observation: [4.07 6.54 1.39 6.17 5.89 4.76 3.67]\n" ] } ], "source": [ "# need to create new columns in order to prepare a lifelines Table\n", "data['observed'] = data['EVENT'] == 1\n", "data['censored'] = data['EVENT'] == 0\n", "\n", "# Transforming survival-table data into lifelines format\n", "from lifelines.utils import survival_events_from_table\n", "\n", "T, E, W = survival_events_from_table(\n", " data.set_index('YEARS'),\n", " observed_deaths_col='observed',\n", " censored_col='censored')\n", "\n", "print(\"Durations of observation:\", T)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Plotting the survival function\n", "\n", "We can then plot the survival function with confidence intervals using `kmf.plot_survival_function()`." ] }, { "cell_type": "code", "execution_count": 50, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkwAAAH2CAYAAACC3dNNAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAABRPklEQVR4nO3deXxOZ/7/8fftjiyyk4UQidpVxVaZaFEVW42itVS1QlW/7YRS1WnNTEWXka6GYmi1llYNXdBOtVTT0GmrRUzU0NpKY1qCFlmUkFy/P/zc01sSJyGcLK/n43EecV9n+5zLfcvbOdd9jsMYYwQAAIBiVbO7AAAAgPKOwAQAAGCBwAQAAGCBwAQAAGCBwAQAAGCBwAQAAGCBwAQAAGCBwAQAAGCBwAQAAGCBwASUEyNGjFB0dPQV3YfD4dCUKVMsl8vMzNTAgQNVq1YtORwOTZ8+/YrWdamuRp9drpycHN17772qXbu2HA6Hxo8fX+ptTJkyRQ6HQ0ePHi37Aiugm266STfddFOZbvN8HwPFITChStq2bZsGDhyoqKgoeXt7q27duurevbtmzpxpd2nlwkMPPaQ1a9Zo0qRJeuONN9SrVy/bavnpp580ZcoUpaen21bD5Zg6daoWLlyoBx54QG+88Ybuvvvuiy67cuXKq1dcFXPy5ElNmTJF69ats7sUVEAOniWHqubLL79U165dVb9+fSUkJKh27do6cOCAvvrqK+3du1d79uyxpa4zZ86ooKBAXl5eV2wfDodDSUlJlmeZateurfj4eC1evPiK1VJSmzdv1vXXX68FCxZoxIgRbvOuRp9drt/97nfy8PDQ559/brmsn5+fBg4cqIULF7q1T5kyRU888YSOHDmikJCQK1RpxXH+7FJpg8/Ro0cVGhpa5Gfg7NmzOnv2rLy9vcumSFQ6HnYXAFxtf/3rXxUYGKhNmzYpKCjIbd7hw4fLbD+5ubny9fUt8fLVq1cvs31frsOHDxfqm/KoPPVZcQ4fPqwWLVrYXQYseHh4yMODX4koHpfkUOXs3btX1157bZGBICwszPXn/fv3y+FwFPrfvlR4LND58Q87duzQnXfeqeDgYN1444164YUX5HA49MMPPxTaxqRJk+Tp6aljx45Jch+Pc+bMGdWsWVMjR44stF5WVpa8vb01ceJESVJeXp4mT56sdu3aKTAwUL6+vurUqZNSU1NL0SvnLFy4UA6HQ8YYzZ49Ww6HwzWuo7gxHufX2b9/v6stOjpav//97/X555+rQ4cO8vb21jXXXKPXX3+90PrHjx/XQw89pOjoaHl5ealevXoaPny4jh49qnXr1un666+XJI0cOdJVz/m/k6LGMOXm5urhhx9WZGSkvLy81LRpU73wwgu68GS6w+HQmDFjtHLlSrVs2VJeXl669tprtXr16hL11eHDhzVq1CiFh4fL29tbMTExWrRokWv+unXr5HA4tG/fPq1atcpV+2/76cJ6cnNztWjRIteyF55RO378uEaMGKGgoCAFBgZq5MiROnnyZKFtLV68WO3atZOPj49q1qypO+64QwcOHLA8pvN/x999950GDx6sgIAA1apVS+PGjdOpU6fclj179qyeeuopNWzYUF5eXoqOjtaf/vQnnT592m258++Fjz/+WK1bt5a3t7datGih5cuXF7nvCxX1/rpQST4D+/fvV2hoqCTpiSeecPXx+c9xUfsv7TGW5P2OiovAhConKipKaWlp+s9//lPm2x40aJBOnjypqVOnavTo0Ro8eLAcDofeeuutQsu+9dZb6tGjh4KDgwvNq169ugYMGKCVK1cqLy/Pbd7KlSt1+vRp3XHHHZLOBahXX31VN910k5599llNmTJFR44cUc+ePUs97qdz58564403JEndu3fXG2+84XpdWnv27NHAgQPVvXt3vfjiiwoODtaIESO0fft21zI5OTnq1KmTZs6cqR49emjGjBm6//779d133+m///2vmjdvrieffFKSdN9997nq6dy5c5H7NMbo1ltv1d/+9jf16tVL06ZNU9OmTfXII49owoQJhZb//PPP9Yc//EF33HGHnnvuOZ06dUq33367fv7554se26+//qqbbrpJb7zxhoYNG6bnn39egYGBGjFihGbMmCFJat68ud544w2FhISodevWrtrP/9K+0BtvvCEvLy916tTJtez//d//uS0zePBgZWdnKzk5WYMHD9bChQv1xBNPuC3z17/+VcOHD1fjxo01bdo0jR8/XikpKercubOOHz9+0eP67X5OnTql5ORk3XLLLXrppZd03333uS1z7733avLkyWrbtq3+9re/qUuXLkpOTna9L39r9+7dGjJkiHr37q3k5GR5eHho0KBBWrt2bYnqsVKSz0BoaKjmzJkjSRowYICrj2+77bZit1uaYyzJ+x0VnAGqmI8//tg4nU7jdDpNXFyc+eMf/2jWrFlj8vLy3Jbbt2+fkWQWLFhQaBuSTFJSkut1UlKSkWSGDh1aaNm4uDjTrl07t7aNGzcaSeb11193tSUkJJioqCjX6zVr1hhJ5p///Kfburfccou55pprXK/Pnj1rTp8+7bbMsWPHTHh4uLnnnnsuWndxJJnExES3tvPHeKEFCxYYSWbfvn2utqioKCPJfPbZZ662w4cPGy8vL/Pwww+72iZPnmwkmeXLlxfabkFBgTHGmE2bNhX793Bhn61cudJIMk8//bTbcgMHDjQOh8Ps2bPH7Rg9PT3d2rZu3WokmZkzZxba129Nnz7dSDKLFy92teXl5Zm4uDjj5+dnsrKy3PqiT58+F93eeb6+viYhIaFQ+/m+v/Dvc8CAAaZWrVqu1/v37zdOp9P89a9/dVtu27ZtxsPDo1B7cfu59dZb3dr/8Ic/GElm69atxhhj0tPTjSRz7733ui03ceJEI8l8+umnrrbz74V3333X1XbixAlTp04d06ZNm0L7vlBR768uXbqYLl26uF6X9DNw5MiRYj8DF+7/Uo7R6v2Oio0zTKhyunfvrg0bNujWW2/V1q1b9dxzz6lnz56qW7eu3n///cva9v3331+obciQIUpLS9PevXtdbcuWLZOXl5f69etX7LZuvvlmhYSEaNmyZa62Y8eOae3atRoyZIirzel0ytPTU5JUUFCgX375RWfPnlX79u21ZcuWyzqey9GiRQt16tTJ9To0NFRNmzbV999/72p79913FRMTowEDBhRa/1K+4v3hhx/K6XTqwQcfdGt/+OGHZYzRRx995NYeHx+vhg0bul63atVKAQEBbjUWt5/atWtr6NChrrbq1avrwQcfVE5OjtavX1/q2kviwvdXp06d9PPPPysrK0uStHz5chUUFGjw4ME6evSoa6pdu7YaN25c4su0iYmJbq/Hjh0r6dxx//bnhWftHn74YUnSqlWr3NojIiLc/o4DAgI0fPhw/fvf/9ahQ4dKVNPFXInPQGmPsSTvd1RsBCZUSddff72WL1+uY8eOaePGjZo0aZKys7M1cOBA7dix45K326BBg0JtgwYNUrVq1VzBxxijt99+W71791ZAQECx2/Lw8NDtt9+u9957zzVmYvny5Tpz5oxbYJKkRYsWqVWrVvL29latWrUUGhqqVatW6cSJE5d8LJerfv36hdqCg4NdY7akc+PJWrZsWWb7/OGHHxQRESF/f3+39ubNm7vml7bG4vbTuHFjVavm/k9ocfspKxfWe/5y7vl6d+/eLWOMGjdurNDQULfp22+/LfGXGho3buz2umHDhqpWrZprHNEPP/ygatWqqVGjRm7L1a5dW0FBQYWOv1GjRoUCcJMmTSTpomOTSqOsPwOlPcZLfS+h4iAwoUrz9PTU9ddfr6lTp2rOnDk6c+aM3n77bUnFn+HIz88vdns+Pj6F2iIiItSpUyfXOKavvvpKGRkZhUJPUe644w5lZ2e7zoy89dZbatasmWJiYlzLLF68WCNGjFDDhg312muvafXq1Vq7dq1uvvlmFRQUWO6jpErbH06ns8h2U47uZFIRavwtq3oLCgrkcDhc74ELp5dffvmS9lvc331Z3ujxUj5v513Jz0BJj7GivZdQenyHEvj/2rdvL0k6ePCgpP/97/3CgbKXcvZgyJAh+sMf/qCdO3dq2bJlqlGjhvr27Wu5XufOnVWnTh0tW7ZMN954oz799FP9+c9/dlvmnXfe0TXXXKPly5e7/eOelJRU6jov5rf98dtvGF7O2ZSGDRtaDr4vzS/lqKgoffLJJ8rOznY7y/Tdd9+55peFqKgoffPNNyooKHA7y3S5+7ncANKwYUMZY9SgQQPXGZxLsXv3brezpXv27FFBQYHrG4lRUVEqKCjQ7t27XWfVpHN3iD9+/Hih49+zZ4+MMW7Ht2vXLklybfNy3l8l/QyU9r1UmmNE5ccZJlQ5qampRf6v7/yYhaZNm0o6N84iJCREn332mdtyf//730u9z9tvv11Op1P/+Mc/9Pbbb+v3v/99ie7RVK1aNQ0cOFD//Oc/9cYbb+js2bOFzkyd/5/tb4/p66+/1oYNG0pd58WcH+vz2/44/zX4S3X77bdr69atWrFiRaF554/nfD+V5Btet9xyi/Lz8zVr1iy39r/97W9yOBzq3bv3Jdd64X4OHTrkNr7s7Nmzmjlzpvz8/NSlS5dL2q6vr2+Jv8lWlNtuu01Op1NPPPFEofe4Mcby23/nzZ492+31+Tvgn++/W265RZIKPTJn2rRpkqQ+ffq4tf/0009uf8dZWVl6/fXX1bp1a9WuXVvS5b2/SvoZqFGjhqSSv5ekkh8jKj/OMKHKGTt2rE6ePKkBAwaoWbNmysvL05dffqlly5YpOjra7d5H9957r5555hnde++9at++vT777DPX/4xLIywsTF27dtW0adOUnZ1dostx5w0ZMkQzZ85UUlKSrrvuOrf/7UrS73//ey1fvlwDBgxQnz59tG/fPs2dO1ctWrRQTk5OqWstTo8ePVS/fn2NGjVKjzzyiJxOp+bPn6/Q0FBlZGRc0jYfeeQRvfPOOxo0aJDuuecetWvXTr/88ovef/99zZ07VzExMWrYsKGCgoI0d+5c+fv7y9fXV7GxsUWOF+vbt6+6du2qP//5z9q/f79iYmL08ccf67333tP48ePdBnhfjvvuu08vv/yyRowYobS0NEVHR+udd97RF198oenTpxcaQ1VS7dq10yeffKJp06YpIiJCDRo0UGxsbInXb9iwoZ5++mlNmjRJ+/fvV//+/eXv7699+/ZpxYoVuu+++1z377qYffv26dZbb1WvXr20YcMGLV68WHfeeafrUnBMTIwSEhL0yiuv6Pjx4+rSpYs2btyoRYsWqX///uratavb9po0aaJRo0Zp06ZNCg8P1/z585WZmakFCxa4lrmc91dJPwM+Pj5q0aKFli1bpiZNmqhmzZpq2bJlkePoSnuMqAJs+GYeYKuPPvrI3HPPPaZZs2bGz8/PeHp6mkaNGpmxY8eazMxMt2VPnjxpRo0aZQIDA42/v78ZPHiwOXz4cLG3FThy5Eix+503b56RZPz9/c2vv/5aaP6FX5E/r6CgwERGRhb5dfnz86dOnWqioqKMl5eXadOmjfnggw+K3N6FdRdHRdxWwBhj0tLSTGxsrPH09DT169c306ZNK/a2AkV9lf7Cr4MbY8zPP/9sxowZY+rWrWs8PT1NvXr1TEJCgjl69Khrmffee8+0aNHCeHh4uN1ioKhjzM7ONg899JCJiIgw1atXN40bNzbPP/+86zYFVscYFRVV5Ff7L5SZmWlGjhxpQkJCjKenp7nuuuuKvPVBaW4r8N1335nOnTsbHx8fI8lVR3Hvr6L63hhj3n33XXPjjTcaX19f4+vra5o1a2YSExPNzp07L7r/8/vZsWOHGThwoPH39zfBwcFmzJgxhd6zZ86cMU888YRp0KCBqV69uomMjDSTJk0yp06dKvL416xZY1q1amW8vLxMs2bNzNtvv11o/yV9f134PirNZ+DLL7807dq1M56enm6fh6Jua1DaY7xQUe93VFw8Sw4AIOnKPLMuOjpaLVu21AcffFAm2wPswhgmAAAACwQmAAAACwQmAAAAC4xhAgAAsMAZJgAAAAsEJgAAAAvcuLIMFRQU6KeffpK/v3+ZPmMJAABcGcYYZWdnKyIiotADtX+LwFSGfvrpJ0VGRtpdBgAAKKUDBw6oXr16xc4nMJWh849DOHDggAICAmyuBgAAWMnKylJkZKTlI40ITGXo/GW4gIAAAhMAABWI1VAaBn0DAABYIDABAABYIDABAABYIDABAABYIDABAABYIDABAABYIDABAABYqJSB6bPPPlPfvn0VEREhh8OhlStXWq6zbt06tW3bVl5eXmrUqJEWLlx4xesEAAAVQ6UMTLm5uYqJidHs2bNLtPy+ffvUp08fde3aVenp6Ro/frzuvfderVmz5gpXCgAAKoJKeafv3r17q3fv3iVefu7cuWrQoIFefPFFSVLz5s31+eef629/+5t69ux5pcq0VlAg/fyzffvH/wQHSx6V8uMCACgBfgNI2rBhg+Lj493aevbsqfHjx9tT0Hk//yyFhdlbA84JDJRmzZLuusvuSgAANqiUl+RK69ChQwoPD3drCw8PV1ZWln799ddi1zt9+rSysrLcJlRSJ05IY8ZIZ8/aXQkAwAYEpsuQnJyswMBA1xQZGWl3SbiSTpyQjh2zuwoAgA0ITJJq166tzMxMt7bMzEwFBATIx8en2PUmTZqkEydOuKYDBw5c6VIBAIANGMMkKS4uTh9++KFb29q1axUXF3fR9by8vOTl5XXlCqtVSzp8+MptHxd39KjUooXdVQAAyoFKGZhycnK0Z88e1+t9+/YpPT1dNWvWVP369TVp0iT9+OOPev311yVJ999/v2bNmqU//vGPuueee/Tpp5/qrbfe0qpVq+w6hHOqVZNCQ+2tAQAAVM5Lcps3b1abNm3Upk0bSdKECRPUpk0bTZ48WZJ08OBBZWRkuJZv0KCBVq1apbVr1yomJkYvvviiXn31VXtvKQAAAMoNhzHG2F1EZZGVlaXAwECdOHFCAQEBdpeDy3XkSOHbOnz/vRQSYk89wJXi4SFdZLwmUJmV9Hd3pbwkB1wxa9dKhGFUNn5+UrduhCbgIghMQGn4+xOYULmcOiXl5HCPMcACgQkoDR8fqUYNu6sAylZent0VAOVepRz0DQAAUJYITAAAABYITAAAABYITAAAABYITAAAABYITAAAABYITAAAABYITAAAABYITAAAABYITAAAABYITAAAABYITAAAABYITAAAABYITAAAABYITAAAABYITAAAABY87C4AAFAOnDxpdwWXzsND8vGxuwpUcgQmAKjKnE4pJ0dKTbW7kkvn5yd160ZowhVFYAKAqszLS4qMlPLz7a7k0pw6dS7wnT1rdyWo5AhMAFDVeXnZXcHlycuzuwJUAQz6BgAAsEBgAgAAsEBgAgAAsEBgAgAAsEBgAgAAsEBgAgAAsEBgAoDy5uefpYULz/0EUC4QmACgvPn5Z2nRIgITUI4QmAAAACwQmAAAACwQmAAAACzwLDmgNE6csLsCVAXZ2f/7efz4ld+fv7/kdF75/QAVGIEJKI0RI+yuAFXJxIlXZz++vtK4cVL37ldnf0AFxCU5AKjqcnOlGTOk/Hy7KwHKLQITUJzgYCkw0O4qgKsjN/d/lwIBFEJgAorj4SHNmkVoAgAwhgm4qLvukvr0kd5++9zAWB8fuytCVbB377nxSy+8IDVseGX2ceIEY/KAUiAwAVY8PKSAgHNTjRp2V4OqwN//fz+DgmwtBcA5XJIDAACwQGACAACwQGACAACwQGACAACwQGACgPKmVi0pIeHcTwDlAt+SA4DyplYtvvIPlDOcYQIAALBAYAIAALBAYAIAALDAGCYAQMV38qTdFVQsHh486qmUKnVgmj17tp5//nkdOnRIMTExmjlzpjp06FDs8tOnT9ecOXOUkZGhkJAQDRw4UMnJyfL29r6KVQMASszplHJypNRUuyupWPz8pG7dCE2lUGkD07JlyzRhwgTNnTtXsbGxmj59unr27KmdO3cqLCys0PJLlizRY489pvnz56tjx47atWuXRowYIYfDoWnTptlwBAAAS15eUmSklJ9vdyUVx6lT50Lm2bN2V1KhVNrANG3aNI0ePVojR46UJM2dO1erVq3S/Pnz9dhjjxVa/ssvv9QNN9ygO++8U5IUHR2toUOH6uuvv76qdQMASsnLy+4KKp68PLsrqHAq5aDvvLw8paWlKT4+3tVWrVo1xcfHa8OGDUWu07FjR6WlpWnjxo2SpO+//14ffvihbrnllmL3c/r0aWVlZblNAACg8qmUZ5iOHj2q/Px8hYeHu7WHh4fru+++K3KdO++8U0ePHtWNN94oY4zOnj2r+++/X3/605+K3U9ycrKeeOKJMq0dAACUP5XyDNOlWLdunaZOnaq///3v2rJli5YvX65Vq1bpqaeeKnadSZMm6cSJE67pwIEDV7FiAABwtVTKM0whISFyOp3KzMx0a8/MzFTt2rWLXOfxxx/X3XffrXvvvVeSdN111yk3N1f33Xef/vznP6tatcLZ0svLS15cOwcAoNKrlGeYPD091a5dO6WkpLjaCgoKlJKSori4uCLXOXnyZKFQ5HQ6JUnGmCtXLAAAKPcq5RkmSZowYYISEhLUvn17dejQQdOnT1dubq7rW3PDhw9X3bp1lZycLEnq27evpk2bpjZt2ig2NlZ79uzR448/rr59+7qCEwAAqJoqbWAaMmSIjhw5osmTJ+vQoUNq3bq1Vq9e7RoInpGR4XZG6S9/+YscDof+8pe/6Mcff1RoaKj69u2rv/71r3YdAgAAKCcchutNZSYrK0uBgYE6ceKEAgIC7C4HZSU7W1q1SgoIkGrUsLsaoGwcPy4NGODetmKFFBRkRzW4mk6elLKypD59JH9/u6uxXUl/d1fKMUwAAABlicAEAABggcAEAABggcAEAABggcAEAABggcAEAABggcAEAABggcAEAABggcAEAABggcAEAABggcAEAABggcAEAABggcAEAABggcAEAABggcAEAABggcAEAABggcAEAABggcAEAABggcAEAABggcAEAABggcAEAABggcAEAABggcAEAABggcAEAABggcAEAABgwcPuAgAA5cSJE3ZXcGn8/SWn0+4qUMkRmAAA54wYYXcFl8bXVxo3Ture3e5KUIkRmAAAFVturjRjhnTzzZxpKo2TJ+2uoHQ8PCQfH/t2b9ueAQD28fc/d2YmN9fuSspGbq6UnS0FBdldSfnndEo5OVJqqt2VlI6fn9Stm22hicAEAFWR03nuMtaMGZUnNKFkvLykyEgpP9/uSkru1KlzIe/sWdtKIDABQFXVvfu5y1jZ2XZXUnonTlTcMVflgZeX3RWUXl6erbsnMAFAVeZ0chkLKAHuwwQAAGCBwAQAAGCBwAQAAGCBwAQAAGCBwAQAAGCBwAQAAGCBwAQAAGCBwAQAAGCBwAQAAGCBwAQAAGCBwAQAAGCBwAQAAGCBwAQAAGCBwAQAAGCBwAQAAGCBwAQAAGCBwAQAAGCBwAQAAGChUgem2bNnKzo6Wt7e3oqNjdXGjRsvuvzx48eVmJioOnXqyMvLS02aNNGHH354laoFAADllYfdBVwpy5Yt04QJEzR37lzFxsZq+vTp6tmzp3bu3KmwsLBCy+fl5al79+4KCwvTO++8o7p16+qHH35QUFDQ1S8eAACUK5U2ME2bNk2jR4/WyJEjJUlz587VqlWrNH/+fD322GOFlp8/f75++eUXffnll6pevbokKTo6+mqWDAAAyqlKeUkuLy9PaWlpio+Pd7VVq1ZN8fHx2rBhQ5HrvP/++4qLi1NiYqLCw8PVsmVLTZ06Vfn5+cXu5/Tp08rKynKbAABA5VMpA9PRo0eVn5+v8PBwt/bw8HAdOnSoyHW+//57vfPOO8rPz9eHH36oxx9/XC+++KKefvrpYveTnJyswMBA1xQZGVmmxwEAAMqHShmYLkVBQYHCwsL0yiuvqF27dhoyZIj+/Oc/a+7cucWuM2nSJJ04ccI1HThw4CpWDAAArpZKOYYpJCRETqdTmZmZbu2ZmZmqXbt2kevUqVNH1atXl9PpdLU1b95chw4dUl5enjw9PQut4+XlJS8vr7ItHgAAlDuV8gyTp6en2rVrp5SUFFdbQUGBUlJSFBcXV+Q6N9xwg/bs2aOCggJX265du1SnTp0iwxIAAKg6KmVgkqQJEyZo3rx5WrRokb799ls98MADys3NdX1rbvjw4Zo0aZJr+QceeEC//PKLxo0bp127dmnVqlWaOnWqEhMT7ToEAABQTlTKS3KSNGTIEB05ckSTJ0/WoUOH1Lp1a61evdo1EDwjI0PVqv0vL0ZGRmrNmjV66KGH1KpVK9WtW1fjxo3To48+atchAACAcqLSBiZJGjNmjMaMGVPkvHXr1hVqi4uL01dffXWFqwIAABVNpb0kBwAAUFYITAAAABYITAAAABYITAAAABYITAAAABYITAAAABYITAAAABYITAAAABYITAAAABYITAAAABYITAAAABYITAAAABYITAAAABYITAAAABYITAAAABYITAAAABYITAAAABY87C4AAIAyceKE3RVUHP7+ktNpdxUVCoEJAFA5jBhhdwUVh6+vNG6c1L273ZVUGFySAwCgqsnNlWbMkPLz7a6kwiAwAQAqHn//c2dJcOlyc6XsbLurqDAITACAisfpPHdJidCEq4QxTACAiql7d+nmmzlLUlInTjDO6zIQmAAAFZfTKQUF2V0FqgAuyQEAAFggMAEAAFggMAEAAFggMAEAAFggMAEAAFggMAEAAFggMAEAAFggMAEAAFggMAEAAFggMAEAAFggMAEAAFggMAEAAFggMAEAAFggMAEAAFggMAEAAFggMAEAAFggMAEAAFggMAEAAFggMAEAAFggMAEAAFggMAEAAFggMAEAAFggMAEAAFggMAEAAFggMAEAAFggMAEAAFggMAEAAFio1IFp9uzZio6Olre3t2JjY7Vx48YSrbd06VI5HA7179//yhYIAAAqhEobmJYtW6YJEyYoKSlJW7ZsUUxMjHr27KnDhw9fdL39+/dr4sSJ6tSp01WqFAAAlHeVNjBNmzZNo0eP1siRI9WiRQvNnTtXNWrU0Pz584tdJz8/X8OGDdMTTzyha6655ipWCwAAyrNKGZjy8vKUlpam+Ph4V1u1atUUHx+vDRs2FLvek08+qbCwMI0aNapE+zl9+rSysrLcJgAAUPlUysB09OhR5efnKzw83K09PDxchw4dKnKdzz//XK+99prmzZtX4v0kJycrMDDQNUVGRl5W3QAAoHyqlIGptLKzs3X33Xdr3rx5CgkJKfF6kyZN0okTJ1zTgQMHrmCVAADALh52F3AlhISEyOl0KjMz0609MzNTtWvXLrT83r17tX//fvXt29fVVlBQIEny8PDQzp071bBhw0LreXl5ycvLq4yrBwAA5U2lPMPk6empdu3aKSUlxdVWUFCglJQUxcXFFVq+WbNm2rZtm9LT013Trbfeqq5duyo9PZ1LbQAAVHGV8gyTJE2YMEEJCQlq3769OnTooOnTpys3N1cjR46UJA0fPlx169ZVcnKyvL291bJlS7f1g4KCJKlQOwAAqHoqbWAaMmSIjhw5osmTJ+vQoUNq3bq1Vq9e7RoInpGRoWrVKuUJNgAAUMYqbWCSpDFjxmjMmDFFzlu3bt1F1124cGHZFwQAACokTrEAAABYIDABAABYIDABAIBL9/PP0sKF535eKb/8Ir3zjlTMzaevBgITAAC4dD//LC1adOUD07vvEpgAAADKMwITAACABQITAACAhUp9HyYAAHARJ05c/jays//38/jxy99eUXJyzv08flwqKJBsuPE0gQkAgKpqxIiy29bEiWW3reLceqt0+LAUGnrl93UBAhMAAKg4vvlGCg4u3F6nzrnpCiEwAQBQFfj7S76+Um6u3ZVcnvj4otuTkqQpU67YbglMAABUBU6nNG6cNGNGxQ5Nn3xS/BmmK4jABJTUqVN2VwBcGU6n5OVldxW4Grp3l26++X8DtcvC3r3nxi+98ILUsGHZbfe3duyQ/vxn6f33pa5dGfQNlEseHpKf37lvaeTl2V0NUPZycqTISEJTVeF0SkFBZbc9f////SzL7f6Wn9+5n0FBtoQlicAEWPPxkbp1k86etbsSoOydPCmlpkr5+XZXApRrBCagJHx87K4AAGAj7vQNAABggcAEAABggcAEAAAuXa1aUkLCuZ9XSs2a0u23S7VrX7l9WGAMEwAAuHS1apXtI1aKUrOmNHCgrYGJM0wAAAAWCEwAAAAWCEwAAAAWCEwAAAAWCEwAAAAWCEwAAAAWCEwAAAAWCEwAAAAWCEwAAAAWCEwAAAAWCEwAAAAWCEwAAAAWCEwAAAAWCEwAAAAWCEwAAAAWCEwAAAAWCEwAAAAWCEwAAAAWCEwAAAAWCEwAAAAWCEwAAAAWCEwAAAAWCEwAAAAWCEwAAAAWCEwAAAAWCEwAAAAWCEwAAAAWCEwAAAAWCEwAAAAWPOwuAABQDpw6ZXcFl87plLy87K4ClVylDkyzZ8/W888/r0OHDikmJkYzZ85Uhw4dilx23rx5ev311/Wf//xHktSuXTtNnTq12OUBoFLw8JD8/KScHCkvz+5qLk1OjhQZSWjCFVVpA9OyZcs0YcIEzZ07V7GxsZo+fbp69uypnTt3KiwsrNDy69at09ChQ9WxY0d5e3vr2WefVY8ePbR9+3bVrVvXhiMAgKvAx0fq1k06e9buSi7NyZNSaqqUn293JajkHMYYY3cRV0JsbKyuv/56zZo1S5JUUFCgyMhIjR07Vo899pjl+vn5+QoODtasWbM0fPjwEu0zKytLgYGBOnHihAICAi6rfgBACWRnS6tWSQEBUo0adleDK+XkSSkrS+rTR/L3L9NNl/R3d6Uc9J2Xl6e0tDTFx8e72qpVq6b4+Hht2LChRNs4efKkzpw5o5o1axa7zOnTp5WVleU2AQCAyqdSBqajR48qPz9f4eHhbu3h4eE6dOhQibbx6KOPKiIiwi10XSg5OVmBgYGuKTIy8rLqBgAA5VOlDEyX65lnntHSpUu1YsUKeXt7F7vcpEmTdOLECdd04MCBq1glAAC4WirloO+QkBA5nU5lZma6tWdmZqp27doXXfeFF17QM888o08++UStWrW66LJeXl7y4lsZAABUepXyDJOnp6fatWunlJQUV1tBQYFSUlIUFxdX7HrPPfecnnrqKa1evVrt27e/GqUCAIAKoFKeYZKkCRMmKCEhQe3bt1eHDh00ffp05ebmauTIkZKk4cOHq27dukpOTpYkPfvss5o8ebKWLFmi6Oho11gnPz8/+fn52XYcAADAfpU2MA0ZMkRHjhzR5MmTdejQIbVu3VqrV692DQTPyMhQtWr/O8E2Z84c5eXlaeDAgW7bSUpK0pQpU65m6QAAoJyptPdhsgP3YQKAq4z7MFUN3IcJAACg/CMwAQAAWCAwAQAAWCAwAQAAWCAwAQAAWCAwAQAAWCAwAQAAWCAwAQAAWCAwAQAAWCAwAQAAWCAwAQAAWCAwAQAAWCAwAQAAWCAwAQAAWCAwAQAAWCAwAQAAWCAwAQAAWPCwuwAAAC7bqVN2V1CxOJ2Sl5fdVVQoBCYAQMXl4SH5+Uk5OVJent3VVBw5OVJkJKGpFAhMAICKy8dH6tZNOnvW7koqjpMnpdRUKT/f7koqFAITAKBi8/GxuwJUAQz6BgAAsEBgAgAAsEBgAgAAsEBgAgAAsEBgAgAAsEBgAgAAsEBgAgAAsEBgAgAAsEBgAgAAsEBgAgAAsEBgAgAAsEBgAgAAsEBgAgAAsEBgAgAAsEBgAgAAsEBgAgAAsEBgAgAAsEBgAgAAsOBhdwEAAMAGp07ZXUHJlYNaCUwAAFQlHh6Sn5+UkyPl5dldTcn5+Z2r3SYEJgAAqhIfH6lbN+nsWbsrKR0Pj3O127V72/YMAADsYWPwqKgY9A0AAGCBwAQAAGCBwAQAAGCBwAQAAGCBwAQAAGCBwAQAAGCBwAQAAGCBwAQAAGCBwAQAAGCBwAQAAGCBR6OUIWOMJCkrK8vmSgAAQEmc/519/nd4cQhMZSg7O1uSFBkZaXMlAACgNLKzsxUYGFjsfIexilQosYKCAv3000/y9/eXw+Eos+1mZWUpMjJSBw4cUEBAQJlttyKjTwqjTwqjT9zRH4XRJ4VVtT4xxig7O1sRERGqVq34kUqcYSpD1apVU7169a7Y9gMCAqrEm7c06JPC6JPC6BN39Edh9ElhValPLnZm6TwGfQMAAFggMAEAAFggMFUAXl5eSkpKkpeXl92llBv0SWH0SWH0iTv6ozD6pDD6pGgM+gYAALDAGSYAAAALBCYAAAALBCYAAAALBCYAAAALBKYKYPbs2YqOjpa3t7diY2O1ceNGu0uyzWeffaa+ffsqIiJCDodDK1eutLskWyUnJ+v666+Xv7+/wsLC1L9/f+3cudPusmw1Z84ctWrVynXTvbi4OH300Ud2l1WuPPPMM3I4HBo/frzdpdhmypQpcjgcblOzZs3sLst2P/74o+666y7VqlVLPj4+uu6667R582a7yyoXCEzl3LJlyzRhwgQlJSVpy5YtiomJUc+ePXX48GG7S7NFbm6uYmJiNHv2bLtLKRfWr1+vxMREffXVV1q7dq3OnDmjHj16KDc31+7SbFOvXj0988wzSktL0+bNm3XzzTerX79+2r59u92llQubNm3Syy+/rFatWtldiu2uvfZaHTx40DV9/vnndpdkq2PHjumGG25Q9erV9dFHH2nHjh168cUXFRwcbHdp5QK3FSjnYmNjdf3112vWrFmSzj2vLjIyUmPHjtVjjz1mc3X2cjgcWrFihfr37293KeXGkSNHFBYWpvXr16tz5852l1Nu1KxZU88//7xGjRpldym2ysnJUdu2bfX3v/9dTz/9tFq3bq3p06fbXZYtpkyZopUrVyo9Pd3uUsqNxx57TF988YX+9a9/2V1KucQZpnIsLy9PaWlpio+Pd7VVq1ZN8fHx2rBhg42Vobw6ceKEpHMBAVJ+fr6WLl2q3NxcxcXF2V2O7RITE9WnTx+3f1Oqst27dysiIkLXXHONhg0bpoyMDLtLstX777+v9u3ba9CgQQoLC1ObNm00b948u8sqNwhM5djRo0eVn5+v8PBwt/bw8HAdOnTIpqpQXhUUFGj8+PG64YYb1LJlS7vLsdW2bdvk5+cnLy8v3X///VqxYoVatGhhd1m2Wrp0qbZs2aLk5GS7SykXYmNjtXDhQq1evVpz5szRvn371KlTJ2VnZ9tdmm2+//57zZkzR40bN9aaNWv0wAMP6MEHH9SiRYvsLq1c8LC7AABlIzExUf/5z3+q/DgMSWratKnS09N14sQJvfPOO0pISND69eurbGg6cOCAxo0bp7Vr18rb29vucsqF3r17u/7cqlUrxcbGKioqSm+99VaVvXRbUFCg9u3ba+rUqZKkNm3a6D//+Y/mzp2rhIQEm6uzH2eYyrGQkBA5nU5lZma6tWdmZqp27do2VYXyaMyYMfrggw+UmpqqevXq2V2O7Tw9PdWoUSO1a9dOycnJiomJ0YwZM+wuyzZpaWk6fPiw2rZtKw8PD3l4eGj9+vV66aWX5OHhofz8fLtLtF1QUJCaNGmiPXv22F2KberUqVPoPxXNmzev8pcqzyMwlWOenp5q166dUlJSXG0FBQVKSUlhPAYkScYYjRkzRitWrNCnn36qBg0a2F1SuVRQUKDTp0/bXYZtunXrpm3btik9Pd01tW/fXsOGDVN6erqcTqfdJdouJydHe/fuVZ06dewuxTY33HBDoduS7Nq1S1FRUTZVVL5wSa6cmzBhghISEtS+fXt16NBB06dPV25urkaOHGl3abbIyclx+x/gvn37lJ6erpo1a6p+/fo2VmaPxMRELVmyRO+99578/f1dY9sCAwPl4+Njc3X2mDRpknr37q369esrOztbS5Ys0bp167RmzRq7S7ONv79/oXFtvr6+qlWrVpUd7zZx4kT17dtXUVFR+umnn5SUlCSn06mhQ4faXZptHnroIXXs2FFTp07V4MGDtXHjRr3yyit65ZVX7C6tfDAo92bOnGnq169vPD09TYcOHcxXX31ld0m2SU1NNZIKTQkJCXaXZoui+kKSWbBggd2l2eaee+4xUVFRxtPT04SGhppu3bqZjz/+2O6yyp0uXbqYcePG2V2GbYYMGWLq1KljPD09Td26dc2QIUPMnj177C7Ldv/85z9Ny5YtjZeXl2nWrJl55ZVX7C6p3OA+TAAAABYYwwQAAGCBwAQAAGCBwAQAAGCBwAQAAGCBwAQAAGCBwAQAAGCBwAQAAGCBwAQAAGCBwASg0li3bp0cDoeOHz9+RfezcOFCBQUFuV5PmTJFrVu3vqL7BGAvAhOACuumm27S+PHjXa87duyogwcPKjAw8KrWMXHiRLeHZAOofHj4LoBKw9PTU7Vr177q+/Xz85Ofn99V3y+Aq4czTAAqpBEjRmj9+vWaMWOGHA6HHA6HFi5c6HZJ7vylsw8++EBNmzZVjRo1NHDgQJ08eVKLFi1SdHS0goOD9eCDDyo/P9+17dOnT2vixImqW7eufH19FRsbq3Xr1hVby4WX5EaMGKH+/fvrhRdeUJ06dVSrVi0lJibqzJkzl7wPAPbiDBOACmnGjBnatWuXWrZsqSeffFKStH379kLLnTx5Ui+99JKWLl2q7Oxs3XbbbRowYICCgoL04Ycf6vvvv9ftt9+uG264QUOGDJEkjRkzRjt27NDSpUsVERGhFStWqFevXtq2bZsaN25covpSU1NVp04dpaamas+ePRoyZIhat26t0aNHl9k+AFw9BCYAFVJgYKA8PT1Vo0YN12W47777rtByZ86c0Zw5c9SwYUNJ0sCBA/XGG28oMzNTfn5+atGihbp27arU1FQNGTJEGRkZWrBggTIyMhQRESHp3Bil1atXa8GCBZo6dWqJ6gsODtasWbPkdDrVrFkz9enTRykpKRo9enSZ7QPA1UNgAlCp1ahRwxWWJCk8PFzR0dFuY47Cw8N1+PBhSdK2bduUn5+vJk2auG3n9OnTqlWrVon3e+2118rpdLpe16lTR9u2bSvTfQC4eghMACq16tWru712OBxFthUUFEiScnJy5HQ6lZaW5hZ4JJVqYPfV2AeAq4fABKDC8vT0dBusXRbatGmj/Px8HT58WJ06dSrTbV/NfQAoW3xLDkCFFR0dra+//lr79+/X0aNHXWdwLkeTJk00bNgwDR8+XMuXL9e+ffu0ceNGJScna9WqVWVQ9dXZB4CyRWACUGFNnDhRTqdTLVq0UGhoqDIyMspkuwsWLNDw4cP18MMPq2nTpurfv782bdqk+vXrl8n2r9Y+AJQdhzHG2F0EAABAecYZJgAAAAsEJgAAAAsEJgAAAAsEJgAAAAsEJgAAAAsEJgAAAAsEJgAAAAsEJgAAAAsEJgAAAAsEJgAAAAsEJgAAAAsEJgAAAAsEJgAAAAsEJgAAAAsEJgAAAAsEJgAAAAsEJgAAAAsEJgAAAAsEJgAAAAsEpv9vxIgR6t+/v91lAACAcqjMAlNRgeOdd96Rt7e3XnzxRY0YMUIOh0P3339/oXUTExPlcDg0YsSIsiqnWPv375fD4VB6erpb+4wZM7Rw4cIrvn+CGQAAFc8VO8P06quvatiwYZozZ44efvhhSVJkZKSWLl2qX3/91bXcqVOntGTJEtWvX/9KlVIigYGBCgoKsrUGAABQPl2RwPTcc89p7NixWrp0qUaOHOlqb9u2rSIjI7V8+XJX2/Lly1W/fn21adOmxNsvKChQcnKyGjRoIB8fH8XExOidd95xzT927JiGDRum0NBQ+fj4qHHjxlqwYIEkqUGDBpKkNm3ayOFw6KabbpJU+MzPTTfdpLFjx2r8+PEKDg5WeHi45s2bp9zcXI0cOVL+/v5q1KiRPvroI9c6+fn5GjVqlKuupk2basaMGa75U6ZM0aJFi/Tee+/J4XDI4XBo3bp1kqQDBw5o8ODBCgoKUs2aNdWvXz/t37+/xH0CAACunDIPTI8++qieeuopffDBBxowYECh+ffcc48rvEjS/Pnz3UJVSSQnJ+v111/X3LlztX37dj300EO66667tH79eknS448/rh07duijjz7St99+qzlz5igkJESStHHjRknSJ598ooMHD7qFtwstWrRIISEh2rhxo8aOHasHHnhAgwYNUseOHbVlyxb16NFDd999t06ePCnpXJCrV6+e3n77be3YsUOTJ0/Wn/70J7311luSpIkTJ2rw4MHq1auXDh48qIMHD6pjx446c+aMevbsKX9/f/3rX//SF198IT8/P/Xq1Ut5eXml6hsAAHAFmDKSkJBgPD09jSSTkpJS5Px+/fqZw4cPGy8vL7N//36zf/9+4+3tbY4cOWL69etnEhISLPdz6tQpU6NGDfPll1+6tY8aNcoMHTrUGGNM3759zciRI4tcf9++fUaS+fe//11kfed16dLF3Hjjja7XZ8+eNb6+vubuu+92tR08eNBIMhs2bCi23sTERHP77bcXux9jjHnjjTdM06ZNTUFBgavt9OnTxsfHx6xZs6bYbQMAgKvDoyzDV6tWrXT06FElJSWpQ4cO8vPzK7RMaGio+vTpo4ULF8oYoz59+rjO/pTEnj17dPLkSXXv3t2tPS8vz3VZ74EHHtDtt9/uOgvUv39/dezY8ZKO5zyn06latWrpuuuuc7WFh4dLkg4fPuxqmz17tubPn6+MjAz9+uuvysvLU+vWrS+6n61bt2rPnj3y9/d3az916pT27t1b6roBAEDZKtPAVLduXb3zzjvq2rWrevXqpY8++qhQCJDOXZYbM2aMpHMBozRycnIkSatWrVLdunXd5nl5eUmSevfurR9++EEffvih1q5dq27duikxMVEvvPBCqfZVvXp1t9cOh8OtzeFwSDp3KU6Sli5dqokTJ+rFF19UXFyc/P399fzzz+vrr7+2PKZ27drpzTffLDQvNDS0VDUDAICyV6aBSZKioqK0fv16V2havXp1odB0fmyOw+FQz549S7X9Fi1ayMvLSxkZGerSpUuxy4WGhiohIUEJCQnq1KmTHnnkEb3wwgvy9PSUdG6Adln74osv1LFjR/3hD39wtV14hsjT07PQvtu2batly5YpLCxMAQEBZV4XAAC4PFfkW3KRkZFat26dDh8+rJ49eyorK8ttvtPp1LfffqsdO3bI6XSWatv+/v6aOHGiHnroIS1atEh79+7Vli1bNHPmTC1atEiSNHnyZL333nvas2ePtm/frg8++EDNmzeXJIWFhcnHx0erV69WZmamTpw4UTYHLalx48bavHmz1qxZo127dunxxx/Xpk2b3JaJjo7WN998o507d+ro0aM6c+aMhg0bppCQEPXr10//+te/tG/fPq1bt04PPvig/vvf/5ZZfQAA4NJcsfsw1atXT+vWrdPRo0eLDE0BAQGXfDblqaee0uOPP67k5GQ1b95cvXr10qpVq1y3DPD09NSkSZPUqlUrde7cWU6nU0uXLpUkeXh46KWXXtLLL7+siIgI9evX7/IO9Df+7//+T7fddpuGDBmi2NhY/fzzz25nmyRp9OjRatq0qdq3b6/Q0FB98cUXqlGjhj777DPVr19ft912m5o3b65Ro0bp1KlTnHECAKAccBhjjN1FAAAAlGc8Sw4AAMBCuQtMGRkZ8vPzK3bKyMiwu0QAAFDFlLtLcmfPnr3oI0Gio6Pl4VHmX+4DAAAoVrkLTAAAAOVNubskZ5ebbrpJ48ePL/NlAQBAxVehAtN///tfeXp6qmXLloXmLVy4UEFBQZe87eXLl+upp566jOouz4gRI+RwOApN1157rW012S06OrrIPklMTLS7NFv9+OOPuuuuu1SrVi35+Pjouuuu0+bNm+0uy1b0ibspU6YU+tw0a9bM7rJsM2fOHLVq1cp1O5u4uDh99NFHdpdlq+TkZF1//fXy9/dXWFiY+vfvr507d9pdVrlWoQLTwoULNXjwYGVlZVk+bqSk8vLyJEk1a9Ys8jEuV8uMGTN08OBB13TgwAHVrFlTgwYNsq0mu23atMmtT9auXStJVbpPjh07phtuuEHVq1fXRx99pB07dujFF19UcHCw3aXZhj4p2rXXXuv2+fn888/tLsk29erV0zPPPKO0tDRt3rxZN998s/r166ft27fbXZpt1q9fr8TERH311Vdau3atzpw5ox49eig3N9fu0sovGx/8WyoFBQXmmmuuMatXrzaPPvqoGT16tGteamqqkeQ2JSUlFbmdpKQkExMTY+bNm2eio6ONw+EwxhjTpUsXM27cONdys2fPNo0aNTJeXl4mLCzM3H777a55Fy77wQcfmICAALN48eIyO94VK1YYh8Nh9u/fX2bbrOjGjRtnGjZsaAoKCuwuxTaPPvqoufHGG+0uo1yhTwo7/+8cihccHGxeffVVu8soNw4fPmwkmfXr19tdSrlVYc4wpaam6uTJk4qPj9ddd92lpUuXupJwx44dNX36dAUEBLj+NzVx4sRit7Vnzx69++67Wr58udLT0wvN37x5sx588EE9+eST2rlzp1avXq3OnTsXua0lS5Zo6NChevPNNzVs2LAyOVZJeu211xQfH6+oqKgy22ZFlpeXp8WLF+uee+5xPfS4Knr//ffVvn17DRo0SGFhYWrTpo3mzZtnd1m2ok+Ktnv3bkVEROiaa67RsGHDuCXL/5efn+/6/REXF2d3OeXG+ceE1axZ0+ZKyjG7E1tJ3XnnnWb8+PGu1zExMWbBggWu1wsWLDCBgYGW20lKSjLVq1c3hw8fdmv/7Vmjd9991wQEBJisrKwit3F+2VmzZpnAwECzbt26Uh/Pxfz444/G6XSaZcuWlel2K7Jly5YZp9NpfvzxR7tLsZWXl5fx8vIykyZNMlu2bDEvv/yy8fb2NgsXLrS7NNvQJ4V9+OGH5q233jJbt241q1evNnFxcaZ+/frF/ptWFXzzzTfG19fXOJ1OExgYaFatWmV3SeVGfn6+6dOnj7nhhhvsLqVcqxCB6dixY8bb29ts3rzZ1fb888+7nYYvTWBq1KhRofbfBqasrCxz3XXXmZCQEHPXXXeZxYsXm9zcXLdl69ata6pXr242btx46QdWjKlTp5patWqZ06dPl/m2K6oePXqY3//+93aXYbvq1aubuLg4t7axY8ea3/3udzZVZD/6xNqxY8dMQEBAlb4Edfr0abN7926zefNm89hjj5mQkBCzfft2u8sqF+6//34TFRVlDhw4YHcp5VqFuCS3ZMkSnTp1SrGxsfLw8JCHh4ceffRRff7559q1a1ept+fr63vR+f7+/tqyZYv+8Y9/qE6dOpo8ebJiYmJ0/Phx1zJt2rRRaGio5s+fL1OGt7Iyxmj+/Pm6++675enpWWbbrch++OEHffLJJ7r33nvtLsV2derUUYsWLdzamjdvXqUvt9An1oKCgtSkSRPt2bPH7lJs4+npqUaNGqldu3ZKTk5WTEyMZsyYYXdZthszZow++OADpaamql69enaXU65ViMD02muv6eGHH1Z6erpr2rp1qzp16qT58+dLOvdhyM/PL7N9enh4KD4+Xs8995y++eYb7d+/X59++qlrfsOGDZWamqr33ntPY8eOLbP9rl+/Xnv27NGoUaPKbJsV3YIFCxQWFqY+ffrYXYrtbrjhhkJf/d21a1eVHutGn1jLycnR3r17VadOHbtLKTcKCgp0+vRpu8uwjTFGY8aM0YoVK/Tpp5+qQYMGdpdU7pX7Z4ykp6dry5YtevPNNwvdR2To0KF68skn9fTTTys6Olo5OTlKSUlRTEyMatSooRo1alzSPj/44AN9//336ty5s4KDg/Xhhx+qoKBATZs2dVuuSZMmSk1N1U033SQPDw9Nnz79Ug/T5bXXXlNsbGyR95qqigoKCrRgwQIlJCTwSBxJDz30kDp27KipU6dq8ODB2rhxo1555RW98sordpdmG/qksIkTJ6pv376KiorSTz/9pKSkJDmdTg0dOtTu0mwxadIk9e7dW/Xr11d2draWLFmidevWac2aNXaXZpvExEQtWbJE7733nvz9/XXo0CFJUmBgoHx8fGyurpyy+ZKgpTFjxpgWLVoUOe/gwYOmWrVq5r333jPGnLsOW6tWrRLdVuBCvx3D9K9//ct06dLFBAcHGx8fH9OqVSu3AdgX3lZgx44dJiwszEyYMOGSjvG848ePGx8fH/PKK69c1nYqkzVr1hhJZufOnXaXUm7885//NC1btjReXl6mWbNmvF8MfXKhIUOGmDp16hhPT09Tt25dM2TIELNnzx67y7LNPffcY6Kiooynp6cJDQ013bp1Mx9//LHdZdlKF9yK5/z02y9TwR3PkgMAALBQIcYwAQAA2InABAAAYIHABAAAYIHABAAAYIHAVAYWLlyooKAgu8sAAABXyBUJTP/+9781aNAghYeHy9vbW40bN9bo0aMv6a7cVck333yjTp06ydvbW5GRkXruuefsLqlcmD17tqKjo+Xt7a3Y2Fht3LjR7pJsR58URp/8z2effaa+ffsqIiJCDodDK1eutLsk29EnhdEnpVPmgemDDz7Q7373O50+fVpvvvmmvv32Wy1evFiBgYF6/PHHy3p3ZSY/P18FBQW27T8rK0s9evRQVFSU0tLS9Pzzz2vKlClV+uZ7krRs2TJNmDBBSUlJ2rJli2JiYtSzZ08dPnzY7tJsQ58URp+4y83NVUxMjGbPnm13KeUGfVIYfVJKZXlTp9zcXBMSEmL69+9f5Pxjx465/rxt2zbTq1cv4+vra8LCwsxdd91ljhw54prfpUsXM3bsWPPII4+Y4OBgEx4e7nYzyoKCApOUlGQiIyONp6enqVOnjhk7dqxr/i+//GLuvvtuExQUZHx8fEyvXr3Mrl27XPPPP6z3vffeM82bNzdOp9Ps27fPnDp1yjz88MMmIiLC1KhRw3To0MGkpqa6HceCBQtMZGSk8fHxMf379zcvvPBCiR78ezF///vfTXBwsNsDdx999FHTtGnTy9puRdehQweTmJjoep2fn28iIiJMcnKyjVXZiz4pjD4pniSzYsUKu8soV+iTwugTa2V6hmnNmjU6evSo/vjHPxY5//w4n+PHj+vmm29WmzZttHnzZq1evVqZmZkaPHiw2/KLFi2Sr6+vvv76az333HN68skntXbtWknSu+++q7/97W96+eWXtXv3bq1cuVLXXXeda90RI0Zo8+bNev/997VhwwYZY3TLLbfozJkzrmVOnjypZ599Vq+++qq2b9+usLAwjRkzRhs2bNDSpUv1zTffaNCgQerVq5d2794tSfr66681atQojRkzRunp6eratauefvrpy+67DRs2qHPnzm4P3O3Zs6d27typY8eOXfb2K6K8vDylpaUpPj7e1VatWjXFx8drw4YNNlZmH/qkMPoEwNVQpg/nOh8qLnzm24VmzZqlNm3aaOrUqa62+fPnKzIyUrt27VKTJk0kSa1atVJSUpIkqXHjxpo1a5ZSUlLUvXt3ZWRkqHbt2oqPj1f16tVVv359dejQwVXH+++/ry+++EIdO3aUJL355puKjIzUypUrNWjQIEnSmTNn9Pe//10xMTGSpIyMDC1YsEAZGRmKiIiQdO6ZTKtXr9aCBQs0depUzZgxQ7169XKFwiZNmujLL7/U6tWrL6vvDh06VOjhh+Hh4a55wcHBl7X9iujo0aPKz8939cN54eHh+u6772yqyl70SWH0CYCroUzPMJkSPmVl69atSk1NlZ+fn2s6H7L27t3rWq5Vq1Zu69WpU8c1JmHQoEH69ddfdc0112j06NFasWKFzp49K0n69ttv5eHhodjYWNe6tWrVUtOmTfXtt9+62jw9Pd32sW3bNuXn56tJkyZuta1fv95V17fffuu2XUmKi4sr0XEDAICKqUzPMJ0/M/Tdd99dNETk5OSob9++evbZZwvNq1OnjuvP1atXd5vncDhcA7MjIyO1c+dOffLJJ1q7dq3+8Ic/6Pnnn9f69etLXK+Pj48cDodbXU6nU2lpaXI6nW7L+vn5lXi7l6J27drKzMx0azv/unbt2ld03+VVSEiInE5nkf1Cn9An59EnAK6GMj3D1KNHD4WEhBT7dfjjx49Lktq2bavt27crOjpajRo1cpt8fX1LvD8fHx/17dtXL730ktatW6cNGzZo27Ztat68uc6ePauvv/7atezPP/+snTt3qkWLFsVur02bNsrPz9fhw4cL1XX+H97mzZu7bVeSvvrqqxLXXJy4uDh99tlnbmOs1q5dq6ZNm1bJy3HSuTOA7dq1U0pKiqutoKBAKSkpVfasHn1SGH0C4Goo08Dk6+urV199VatWrdKtt96qTz75RPv379fmzZv1xz/+Uffff78kKTExUb/88ouGDh2qTZs2ae/evVqzZo1Gjhyp/Pz8Eu1r4cKFeu211/Sf//xH33//vRYvXiwfHx9FRUWpcePG6tevn0aPHq3PP/9cW7du1V133aW6deuqX79+xW6zSZMmGjZsmIYPH67ly5dr37592rhxo5KTk7Vq1SpJ0oMPPqjVq1frhRde0O7duzVr1qzLHr8kSXfeeac8PT01atQobd++XcuWLdOMGTM0YcKEy952RTZhwgTNmzdPixYt0rfffqsHHnhAubm5GjlypN2l2YY+KYw+cZeTk6P09HSlp6dLkvbt26f09HRlZGTYW5iN6JPC6JNSuhJfvdu0aZO57bbbTGhoqPHy8jKNGjUy9913n9m9e7drmV27dpkBAwa4vvbfrFkzM378eFNQUGCMOXdbgXHjxrltt1+/fiYhIcEYY8yKFStMbGysCQgIML6+vuZ3v/ud+eSTT1zLnr+tQGBgoPHx8TE9e/Ys8rYCF8rLyzOTJ0820dHRpnr16qZOnTpmwIAB5ptvvnEt89prr5l69eoZHx8f07dv3zK5rYAxxmzdutXceOONxsvLy9StW9c888wzl73NymDmzJmmfv36xtPT03To0MF89dVXdpdkO/qkMPrkf1JTU42kQtP5fz+rIvqkMPqkdBzGlHCkNgAAQBXFs+QAAAAsEJgAAAAsEJgAAAAsEJgAAAAsEJgAAAAslNvAtH//fjkcjiKnsrhRZGk4HA6tXLnyiu8nIyNDffr0UY0aNRQWFqZHHnnE9biXqmz27NmKjo6Wt7e3YmNjtXHjRrtLss1nn32mvn37KiIi4qq9L8s7+qQw+sTdnDlz1KpVKwUEBCggIEBxcXH66KOP7C6rXHnmmWfkcDg0fvx4u0spt8ptYDrvk08+0cGDB92mdu3a2V1WmcvPz1efPn2Ul5enL7/8UosWLdLChQs1efJku0uz1bJlyzRhwgQlJSVpy5YtiomJUc+ePV3PFKxqcnNzFRMTo9mzZ9tdSrlBnxRGn7irV6+ennnmGaWlpWnz5s26+eab1a9fP23fvt3u0sqFTZs26eWXXy70/FZcwO4bQRVn3759RpL597//XeT8nTt3Gknm22+/dWufNm2aueaaa1yvt23bZnr16mV8fX1NWFiYueuuu8yRI0dc87t06WLGjh1rHnnkERMcHGzCw8NNUlKSa35UVJTbDb2ioqKMMcakp6ebm266yfj5+Rl/f3/Ttm1bs2nTpks+3g8//NBUq1bNHDp0yNU2Z84cExAQYE6fPn3J263oOnToYBITE12v8/PzTUREhElOTraxqvJBklmxYoXdZZQr9Elh9EnRgoODzauvvmp3GbbLzs42jRs3NmvXri3yhtH4n3J/hqk4TZo0Ufv27fXmm2+6tb/55pu68847JZ17dt3NN9+sNm3aaPPmzVq9erUyMzM1ePBgt3UWLVokX19fff3113ruuef05JNPau3atZLOJW9JWrBggQ4ePOh6PWzYMNWrV0+bNm1SWlqaHnvssUIPCy6NDRs26LrrrlN4eLirrWfPnsrKyqqy/wvKy8tTWlqa4uPjXW3VqlVTfHy8NmzYYGNlACqq/Px8LV26VLm5uTxrUOceVdanTx+3f2dRNA+7C7DSsWNHVavmnutycnIknQsts2bN0lNPPSVJ2rVrl9LS0rR48WJJ0qxZs9SmTRtNnTrVte78+fMVGRmpXbt2qUmTJpKkVq1aKSkpSZLUuHFjzZo1SykpKerevbtCQ0MlSUFBQW5PPs/IyNAjjzyiZs2auda7HIcOHXILS5Jcrw8dOnRZ266ojh49qvz8/CL75bvvvrOpKgAV0bZt2xQXF6dTp07Jz89PK1asuOjD2KuCpUuXasuWLa4TAbi4cn+GadmyZa6HA/72IYGSdMcdd2j//v2uQeBvvvmm2rZt6woxW7duVWpqqvz8/FzT+Xl79+51befC67Z16tSxHCMzYcIE3XvvvYqPj9czzzzjtj0AQPnStGlTpaen6+uvv9YDDzyghIQE7dixw+6ybHPgwAGNGzdOb775pry9ve0up0Io94EpMjJSjRo1cpvOq127tm6++WYtWbJEkrRkyRINGzbMNT8nJ0d9+/YtFLh2796tzp07u5a78FKaw+FQQUHBReuaMmWKtm/frj59+ujTTz9VixYttGLFiks+ztq1ayszM9Ot7fzr357ZqkpCQkLkdDqL7Jeq2icALo2np6caNWqkdu3aKTk5WTExMZoxY4bdZdkmLS1Nhw8fVtu2beXh4SEPDw+tX79eL730kjw8PJSfn293ieVOuQ9MVoYNG6Zly5Zpw4YN+v7773XHHXe45rVt21bbt29XdHR0odDl6+tb4n1Ur169yDdPkyZN9NBDD+njjz/WbbfdpgULFlzyccTFxWnbtm1uZ7bWrl2rgICAKnva2NPTU+3atVNKSoqrraCgQCkpKYw9AHBZCgoKdPr0abvLsE23bt20bds2t5MJ7du317Bhw5Seni6n02l3ieVOuR/D9PPPPxcawxMUFOQ6hXjbbbfpgQce0AMPPKCuXbsqIiLCtVxiYqLmzZunoUOH6o9//KNq1qypPXv2aOnSpXr11VdL/IaIjo5WSkqKbrjhBnl5ecnb21uPPPKIBg4cqAYNGui///2vNm3apNtvv/2Sj7NHjx5q0aKF7r77bj333HM6dOiQ/vKXvygxMVFeXl6XvN2KbsKECUpISFD79u3VoUMHTZ8+Xbm5uRo5cqTdpdkiJydHe/bscb3et2+f0tPTVbNmTdWvX9/GyuxDnxRGn7ibNGmSevfurfr16ys7O1tLlizRunXrtGbNGrtLs42/v79atmzp1ubr66tatWoVasf/Z/fX9Ipz/rYCRU3/+Mc/3JYdPHiwkWTmz59faDu7du0yAwYMMEFBQcbHx8c0a9bMjB8/3hQUFBhjTJFfo+zXr59JSEhwvX7//fdNo0aNjIeHh4mKijKnT582d9xxh4mMjDSenp4mIiLCjBkzxvz666+Xdcz79+83vXv3Nj4+PiYkJMQ8/PDD5syZM5e1zcpg5syZpn79+sbT09N06NDBfPXVV3aXZJvU1NQiPxO/fb9WNfRJYfSJu3vuucdERUUZT09PExoaarp162Y+/vhju8sqd7itwMU5jDHmagY0AACAiqbCj2ECAAC40ghMAAAAFghMAAAAFghMAAAAFghMAAAAFghMAAAAFghMAAAAFghMAAAAFghMAAAAFghMAAAAFghMAAAAFghMAAAAFghMAAAAFghMAAAAFghMAAAAFghMAAAAFghMAAAAFghMAAAAFghMAAAAFghMAAAAFghMAAAAFghMAAAAFghMAAAAFghMAAAAFghMAAAAFghMAAAAFghMAAAAFghMAAAAFv4f70L5hIWDOM8AAAAASUVORK5CYII=", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# let's see how it looks\n", "plt.figure(figsize=(6,5))\n", "\n", "kmf.plot(\n", " show_censors=True,\n", " ci_show=True,\n", " legend=False,\n", " lw=3,\n", " c='red',\n", " at_risk_counts=True, # adds summary tables under the plot\n", ")\n", "plt.ylabel('Suvival probability estimate')\n", "plt.title('Survival function of the population');" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "`kmf.plot_survival_function()` plots the Kaplan-Meier estimate of the survival function. The y-axis represents the probability that the event of interest has not yet occurred, while the x-axis represents time.\n", "\n", "The shaded area around the survival function line represents the confidence interval. By default, it's a 95% confidence interval, meaning we are 95% confident that the true survival function lies within this area." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Kaplan-Meier estimator\n", "\n", "### Fraction surviving\n", "\n", "The **fraction surviving**, often represented as S(t), is a fundamental concept in survival analysis. It represents the probability that an individual survives from the time origin (t=0) to a specified future time t.\n", "\n", "The survival function S(t), also known as the probability that life is longer than t, is given by $S(t) = P(T > t)$, where T is the random lifetime taken from the population.\n", "\n", "In the context of survival analysis, \"surviving\" doesn't necessarily mean \"not dying.\" Instead, it means \"not experiencing the event of interest.\" The event could be anything: death, failure of a machine, recovery from a disease, etc.\n", "\n", "Now, let's consider the example of the 7 individuals in the previous study. Let's focus on the 3 first events:\n", "\n", "1. At time t1, one event (death) occurs. So, the fraction surviving drops by 1/7 (since all 7 individuals were at risk at this time). The survival probability at time t1 is thus $S(t_1) = 6/7$.\n", "2. At time t2, one individual is censored. This means they leave the study or the study ends before they have the event. Censoring does not affect the survival probability, so $S(t_2) = S(t_1) = 6/7$.\n", "3. At time t3, another event occurs. Now, only 5 individuals were at risk (since one event occurred at t1, reducing the at-risk group to 6, and one individual was censored at t2, further reducing the at-risk group to 5), so the fraction surviving drops by 1/5. The survival fraction at time t3 is thus $S(t_3) = S(t_2) \\times (1 - 1/5) = (6/7)(4/5)$. And so on.\n", "\n", "We can compute the surviving fractions manually as follows:" ] }, { "cell_type": "code", "execution_count": 51, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0.8571428571428571\n", "0.6857142857142857\n", "0.5142857142857142\n" ] } ], "source": [ "print(6/7)\n", "print(6/7 * 4/5)\n", "print(6/7 * 4/5 * 3/4)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Survival function\n", "\n", "In survival analysis, the survival function S(t) is often estimated using the product of conditional survival probabilities. This is the basis for the Kaplan-Meier estimator, which is one of the most commonly used methods in survival analysis.\n", "\n", "The Kaplan-Meier estimate of the survival function at time t is given by:\n", "\n", "$$\\hat{S}(t) = \\prod_{k \\mid t_k \\le{t}} \\left(1 - \\frac{d_k}{n_k} \\right)$$\n", "\n", "where:\n", "\n", "- The product is over all times $t_i$ up to and including $t$ where at least one event occurred.\n", "- $d_k$ is the number of events that occurred at time $t_k$.\n", "- $n_k$ is the number of individuals at risk at time $t_k$, i.e., the number of individuals who have not yet had an event and have not been censored before time $t_k$.\n", "\n", "This product ensures that the survival function estimate is properly adjusted at each time an event occurs, taking into account the decreasing number of individuals at risk.\n", "\n", "Also note that\n", "\n", "$$\\hat{S}(t_j) = \\prod_{i=1}^j\\frac{n_i - d_i}{n_i} = \\frac{n_j - d_j}{n_j} \\times \\prod_{i=1}^{j-1} \\frac{n_i - d_i}{n_i} = \\frac{n_j - d_j}{n_j} \\times \\hat{S}(t_{j-1})$$" ] }, { "cell_type": "code", "execution_count": 52, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Time: 1.39, Survival fraction S(t): 0.8571428571428572\n", "Time: 3.67, Survival fraction S(t): 0.8571428571428572\n", "Time: 4.07, Survival fraction S(t): 0.6857142857142858\n", "Time: 4.76, Survival fraction S(t): 0.5142857142857143\n", "Time: 5.89, Survival fraction S(t): 0.3428571428571429\n", "Time: 6.17, Survival fraction S(t): 0.3428571428571429\n", "Time: 6.54, Survival fraction S(t): 0.3428571428571429\n" ] } ], "source": [ "import numpy as np\n", "\n", "# convert back from Series to array\n", "durations = data['YEARS'].to_numpy()\n", "event_observed = data['EVENT'].to_numpy()\n", "\n", "# Sort durations and event_observed by durations\n", "sort_indices = np.argsort(durations)\n", "durations_sorted = durations[sort_indices]\n", "events_sorted = event_observed[sort_indices]\n", "\n", "n_at_risk = len(durations_sorted) # initial number at risk\n", "survival_prob = 1.0 # initial survival probability\n", "\n", "for i in range(len(durations_sorted)):\n", " if events_sorted[i]: # if event occurred\n", " survival_prob *= (1 - 1/n_at_risk)\n", " print(f\"Time: {durations_sorted[i]}, Survival fraction S(t): {survival_prob}\")\n", " n_at_risk -= 1 # decrease number at risk" ] }, { "cell_type": "code", "execution_count": 53, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Time: 4.07, Survival fraction S(t): 0.6857142857142858\n" ] } ], "source": [ "# for example for t=4.07, since no event were observed\n", "# at the second time point t=3.67 (censored)\n", "S3 = (1 - 1/7) * (1 - 0/6) * (1 - 1/5)\n", "\n", "print(f\"Time: {4.07}, Survival fraction S(t): {S3}\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We won't rewrite the entire `lifelines` library, but here is an example of how we could have written the Kaplan-Meier survival function:" ] }, { "cell_type": "code", "execution_count": 54, "metadata": {}, "outputs": [], "source": [ "def kaplan_meier(durations, event_observed):\n", " \"\"\"\n", " Estimates the survival function using the Kaplan-Meier method.\n", "\n", " Args:\n", " durations: A numpy array of event times.\n", " event_observed: A numpy array of indicators (1 for event, 0 for censoring) \n", " at each time point.\n", "\n", " Returns:\n", " A tuple containing three numpy arrays:\n", " - survival_probsorted: Estimated survival probabilities at each time point (sorted).\n", " - events_sorted: Unique event (sorted).\n", " - durations_sorted: Unique event times (sorted).\n", " \n", " \"\"\"\n", "\n", " # Sort durations and event_observed by durations\n", " sort_indices = np.argsort(durations)\n", " durations_sorted = np.array(durations)[sort_indices]\n", " events_sorted = np.array(event_observed)[sort_indices]\n", "\n", " # Initialize variables\n", " n_at_risk = len(durations_sorted) # initial number at risk\n", " survival_prob_sorted = np.ones_like(durations_sorted) # initial survival probabilities\n", "\n", " for i in range(len(durations_sorted)):\n", " if events_sorted[i] == 1: # if event occurred\n", " survival_prob_sorted[i:] *= (1 - 1/n_at_risk)\n", " n_at_risk -= 1 # decrease number at risk\n", "\n", " return survival_prob_sorted, events_sorted, durations_sorted" ] }, { "cell_type": "code", "execution_count": 55, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "(array([0.85714286, 0.85714286, 0.68571429, 0.51428571, 0.34285714,\n", " 0.34285714, 0.34285714]), array([1, 0, 1, 1, 1, 0, 0], dtype=int64), array([1.39, 3.67, 4.07, 4.76, 5.89, 6.17, 6.54]))\n" ] } ], "source": [ "print(kaplan_meier(data['YEARS'], data['EVENT']))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Confidence interval\n", "\n", "### The basics\n", "\n", "The confidence interval for the survival function at each time point is typically calculated in using the [Greenwood's formula](https://bookdown.org/sestelo/sa_financial/pointwise-confidence-interval-for-st.html), which provides an estimate of the variance of the survival function. The formula for the variance at time $t$ using the so-called \"linear\" or \"plain\" method is:\n", "\n", "$$s^2_{S(t)} = S^2(t) \\times \\sum_{t_k \\le{t}} \\frac{d_k}{n_k (n_k - d_k)}$$\n", "\n", "where:\n", "\n", "- $S(t)$ is the Kaplan-Meier estimate of the survival function up to time $t$.\n", "- The sum is over all times $t_k$ up to and including $t$ where at least one event occurred.\n", "- $d_k$ is the number of events that occurred at time $t_k$\n", "- $n_k$ is the number of individuals at risk at time $t_k$.\n", "\n", "Therefore, assuming that the sample follow a Normal distribution, so that $S(t) \\pm z^\\ast \\times s_{S(t)}$, the confidence interval, can be determinted:\n", "\n", "$$\\text{IC}^{0.95}_{S(t)} = S(t) \\pm 1.96 \\times s_{S(t)}$$\n", "\n", "where $z^\\ast$ is the z-score or critical value for the desired confidence level (for example, z = 1.96 for a 95% confidence interval), and $s_{S(t)}$ the *standard error* (it will be discussed in more details in later chapters). The common practice is to clip the interval at $[0,1]$.\n", "\n", "Note that these calculations assume that the number of events at each time point and the number of individuals at risk are _large enough_ for the Central Limit Theorem to hold. If these numbers are small, the confidence intervals may not be accurate." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Step-by-step example\n", "\n", "The z-values, also known as **z-scores**, correspond to the number of standard deviations away from the mean in a standard normal distribution. They are used to calculate confidence intervals in statistics, notably using **critical z-values**. These values are derived from the standard normal distribution and are used to calculate the margin of error." ] }, { "cell_type": "code", "execution_count": 56, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Z-value for 99% confidence = 2.575829\n", "Z-value for 95% confidence = 1.959964\n" ] } ], "source": [ "from scipy.stats import norm\n", "\n", "# For a 99% confidence level\n", "z_99 = norm.ppf(0.995) # Two-tailed test, so we use 0.995 instead of 0.99\n", "print(f\"Z-value for 99% confidence = {z_99:.6f}\")\n", "\n", "# For a 95% confidence level\n", "z_95 = norm.ppf(0.975) # Two-tailed test, so we use 0.975 instead of 0.95\n", "print(f\"Z-value for 95% confidence = {z_95:.6f}\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's take the first events in our example to show how we construct the margin of error using the Greenwood's formula. We simply enter the values for $n$ and $d$." ] }, { "cell_type": "code", "execution_count": 57, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "S(t1) = 0.8571\n", "Greenwood variance component at t1 = 0.024\n", "Cumulative sum up to t1 = 0.024\n", "s_S(t1) = 0.132\n", "The confidence interval (plain type) limits at t1 are [0.598 - 1.116]\n" ] } ], "source": [ "# breakdown calculation for the first event\n", "S1 = (1 - 1/7)\n", "component_t1 = 1/(7*(7-1))\n", "cumsum_component_t1 = component_t1\n", "SE1 = S1 * np.sqrt(cumsum_component_t1)\n", "print(f\"S(t1) = {S1:.4f}\")\n", "print(f\"Greenwood variance component at t1 = {component_t1:.3f}\")\n", "print(f\"Cumulative sum up to t1 = {cumsum_component_t1:.3f}\")\n", "print(f\"s_S(t1) = {SE1:.3f}\")\n", "print(f\"The confidence interval (plain type) limits at t1 are \\\n", "[{S1 - 1.96*SE1:.3f} - {S1 + 1.96*SE1:.3f}]\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "When no event is observed, $d=0$ and $e=0$, so that the sum of $e$ doesn't change nor the survival function or its standard error." ] }, { "cell_type": "code", "execution_count": 58, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "S(t3) = 0.6857\n", "Greenwood variance component at t3 = 0.050\n", "Cumulative sum up to t3 = 0.074\n", "s_S(t3) = 0.186\n", "The confidence interval (plain type) limits at t3 are [0.321 - 1.051]\n" ] } ], "source": [ "# breakdown calculation for the third event\n", "S3 = (1 - 1/7) * (1 - 0/6) * (1 - 1/5)\n", "component_t3 = 1 / (5*(5 - 1))\n", "cumsum_component_t3 = component_t1 + 0 + component_t3\n", "SE3 = S3 * np.sqrt(cumsum_component_t3)\n", "print(f\"S(t3) = {S3:.4f}\")\n", "print(f\"Greenwood variance component at t3 = {component_t3:.3f}\")\n", "print(f\"Cumulative sum up to t3 = {cumsum_component_t3:.3f}\")\n", "print(f\"s_S(t3) = {SE3:.3f}\")\n", "print(f\"The confidence interval (plain type) limits at t3 are \\\n", "[{S3 - 1.96*SE3:.3f} - {S3 + 1.96*SE3:.3f}]\")" ] }, { "cell_type": "code", "execution_count": 59, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "S(t4) = 0.5143\n", "Greenwood variance component at t4 = 0.083\n", "Cumulative sum up to t4 = 0.157\n", "s_S(t4) = 0.204\n", "The confidence interval (plain type) limits at t4 are [0.115 - 0.914]\n" ] } ], "source": [ "# breakdown calculation for the 4th event\n", "S4 = (1 - 1/7) * (1 - 0/6) * (1 - 1/5) * (1 - 1/4)\n", "component_t4 = 1 / (4* (4 - 1))\n", "cumsum_component_t4 = component_t1 + 0 + component_t3 + component_t4\n", "SE4 = S4 * np.sqrt(cumsum_component_t4)\n", "print(f\"S(t4) = {S4:.4f}\")\n", "print(f\"Greenwood variance component at t4 = {component_t4:.3f}\")\n", "print(f\"Cumulative sum up to t4 = {cumsum_component_t4:.3f}\")\n", "print(f\"s_S(t4) = {SE4:.3f}\")\n", "print(f\"The confidence interval (plain type) limits at t4 are \\\n", "[{S4 - 1.96*SE4:.3f} - {S4 + 1.96*SE4:.3f}]\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "So that we obtain the following table:\n", "\n", "| time to event | at risk (n) | observed (d) | S(t) | d/(n(n-d)) | cumsum | s_S(t) | LCL | UCL |\n", "| ------------- | ----------- | ------------ | ------ | ---------- | ------ | ------ | ----- | ----- |\n", "| 1.39 | 7 | 1 | 0.8571 | 0.024 | 0.024 | 0.132 | 0.598 | 1.116 |\n", "| 3.67 | 6 | 0 | 0.8571 | 0 | 0.024 | 0.132 | 0.598 | 1.116 |\n", "| 4.07 | 5 | 1 | 0.6857 | 0.050 | 0.074 | 0.186 | 0.321 | 1.051 |\n", "| 4.76 | 4 | 1 | 0.5143 | 0.083 | 0.157 | 0.204 | 0.115 | 0.914 |" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Log survival\n", "\n", "In survival analysis, confidence intervals are often calculated on the log scale, where the _complementary log-log transformation_ $v(t) = \\log(-\\log S(t))$ is used a few reasons:\n", "\n", "1. Symmetry: survival probabilities are bound between 0 and 1, and their distribution can be skewed, especially when probabilities are close to the boundaries. Taking the logarithm helps to symmetrize the distribution, which is a key assumption for the calculation of confidence intervals.\n", "2. Stabilizing variance: the variance of survival probabilities can change over time, and it can be particularly high when probabilities are close to 0 or 1. The log transformation can help to stabilize the variance, making the statistical analysis more reliable.\n", "3. Multiplicative effects: in many cases, particularly in the context of survival analysis, the effects are multiplicative rather than additive. The log transformation converts *multiplicative effects* on the original scale to *additive effects* on the log scale, which simplifies the analysis: $\\log {S(t)} = \\log \\prod (1 - \\frac{d_k}{n_k}) = \\sum{\\log (1 - \\frac{d_k}{n_k})}$.\n", "4. Avoiding impossible values: when calculating confidence intervals for survival probabilities, it's possible to get values that fall outside the range $[0, 1]$, which doesn't make sense for probabilities. The log transformation followed by exponentiation (when calculating the confidence interval bounds) ensures that the confidence intervals fall within the appropriate range.\n", "\n", "The variance for the transformed values is calculated using a specialized formula derived from Greenwood's formula and a technique called the [delta method](https://en.wikipedia.org/wiki/Delta_method). While the derivation involves some calculus (function derivatives), the resulting formula is straightforward to use:\n", "\n", "$$\n", "s^2_{v(t)} = \\frac{\\sum_{t_k \\le t} \\frac{d_k}{n_k(n_k - d_k)}}{\\left( \\sum_{t_k \\le t} \\log \\frac{n_k - d_k}{n_k} \\right)^2}\n", "$$\n", "\n", "where $t_k$ are the distinct event times, $d_k$ is the number of events at time $t_k$, and $n_k$ is the number of individuals at risk at time $t_k$. Note that the numerator in the formula for $s_{v(t)}$ is identical to the summation term in Greenwood's formula, which represents the cumulative sum of estimated hazard increments. So that the confidence interval on the log-log scale is $\\text{CI}_{v(t)} = v(t) \\pm z^\\ast \\times s_{v(t)}$.\n", "\n", "Finally, we can transform back to the original scale, i.e., the scale of the survival function, by exponentiating the lower and upper limits of the confidence interval for the log survival function:\n", "\n", "$$\n", "\\begin{aligned}\n", "\\text{IC}^{0.95}_{S(t)} &= \\exp (-\\exp \\left(\\text{CI}^{0.95}_{v(t)} \\right)) \\\\\n", "&= \\exp (-\\exp \\left( log(-log(S(t))) \\pm 1.96 \\times s_{v(t)} \\right)) \\\\\n", "&= S(t)^{\\exp(\\pm 1.96 \\times s_{v(t)})}\n", "\\end{aligned}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's fill the formula with the values from the first events in our example." ] }, { "cell_type": "code", "execution_count": 60, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Log-survival increment at t1 = -0.154\n", "Log-log variance scaling factor at t1 = -0.154\n", "s_v(t1) = 1.001\n", "The confidence interval (log-log transformation) limits at t1 are [0.3340 - 0.9786]\n" ] } ], "source": [ "# breakdown calculation for the first and second events\n", "denominator_component_t1 = np.log((7-1)/7)\n", "cumsum_denominator_t1 = denominator_component_t1\n", "s_vt1 = np.sqrt(cumsum_component_t1 / cumsum_denominator_t1**2)\n", "print(f\"Log-survival increment at t1 = {denominator_component_t1:.3f}\")\n", "print(f\"Log-log variance scaling factor at t1 = {cumsum_denominator_t1:.3f}\")\n", "print(f\"s_v(t1) = {s_vt1:.3f}\")\n", "print(f\"The confidence interval (log-log transformation) limits at t1 are \\\n", "[{S1**np.exp(1.96*s_vt1):.4f} - {S1**np.exp(-1.96*s_vt1):.4f}]\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Note that when $d=0$, then $l=\\log(n/n)=\\log(1)=0$, so the all the values for the second time-point are identical to the first ones." ] }, { "cell_type": "code", "execution_count": 61, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Log-survival increment at t3 = -0.223\n", "Log-log variance scaling factor at t3 = -0.377\n", "s_v(t3) = 0.720\n", "The confidence interval (log-log transformation) limits at t3 are [0.2128 - 0.9121]\n" ] } ], "source": [ "# breakdown calculation for the third event\n", "denominator_component_t3 = np.log((5-1)/5)\n", "cumsum_denominator_t3 = denominator_component_t1 + 0 + denominator_component_t3\n", "s_vt3 = np.sqrt(cumsum_component_t3 / cumsum_denominator_t3**2)\n", "print(f\"Log-survival increment at t3 = {denominator_component_t3:.3f}\")\n", "print(f\"Log-log variance scaling factor at t3 = {cumsum_denominator_t3:.3f}\")\n", "print(f\"s_v(t3) = {s_vt3:.3f}\")\n", "print(f\"The confidence interval (log-log transformation) limits at t3 are \\\n", "[{S3**np.exp(1.96*s_vt3):.4f} - {S3**np.exp(-1.96*s_vt3):.4f}]\")" ] }, { "cell_type": "code", "execution_count": 62, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Log-survival increment at t4 = -0.288\n", "Log-log variance scaling factor at t4 = -0.665\n", "s_v(t4) = 0.596\n", "The confidence interval (log-log transformation) limits at t4 are [0.1178 - 0.8133]\n" ] } ], "source": [ "# breakdown calculation for the 4th event\n", "denominator_component_t4 = np.log((4-1)/4)\n", "cumsum_denominator_t4 = denominator_component_t1 + 0 + denominator_component_t3 + denominator_component_t4\n", "s_vt4 = np.sqrt(cumsum_component_t4 / cumsum_denominator_t4**2)\n", "print(f\"Log-survival increment at t4 = {denominator_component_t4:.3f}\")\n", "print(f\"Log-log variance scaling factor at t4 = {cumsum_denominator_t4:.3f}\")\n", "print(f\"s_v(t4) = {s_vt4:.3f}\")\n", "print(f\"The confidence interval (log-log transformation) limits at t4 are \\\n", "[{S4**np.exp(1.96*s_vt4):.4f} - {S4**np.exp(-1.96*s_vt4):.4f}]\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "These calculations lead to the final table:\n", "\n", "| time to event | n | d | S(t) | e | cumsum(e) | log((n-d)/n) | cumsum | s_v(t) | LCL | UCL |\n", "| ------------- | - | - | ------ | ----- | --------- | ------------ | ------ | ------ | ------ | ------ |\n", "| 1.39 | 7 | 1 | 0.8571 | 0.024 | 0.024 | -0.154 | -0.154 | 1.001 | 0.3340 | 0.9786 |\n", "| 3.67 | 6 | 0 | 0.8571 | 0 | 0.024 | 0 | -0.154 | 1.001 | 0.3340 | 0.9786 |\n", "| 4.07 | 5 | 1 | 0.6857 | 0.050 | 0.074 | -0.223 | -0.377 | 0.720 | 0.2128 | 0.9121 |\n", "| 4.76 | 4 | 1 | 0.5143 | 0.083 | 0.157 | -0.288 | -0.665 | 0.596 | 0.1178 | 0.8133 |" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Process automation with scripts\n", "\n", "While libraries like `lifelines` offer convenient functions for various statistical analyses, understanding the underlying calculations is valuable. This section dives into creating one script to calculate confidence intervals." ] }, { "cell_type": "code", "execution_count": 63, "metadata": {}, "outputs": [], "source": [ "from scipy.stats import norm\n", "\n", "def greenwood_confidence_interval(\n", " survival_prob_sorted, event_indicators_sorted, alpha=.05, method=\"log-log\"):\n", " \"\"\"\n", " Calculates the confidence intervals, plain and log-log, for a Kaplan-Meier curve, \n", " assuming pre-sorted data, using the Greenwood variance estimate.\n", "\n", " Args:\n", " survival_prob_sorted: A numpy array of estimated survival probabilities \n", " at each time point (assumed to be sorted).\n", " event_indicators: A numpy array of indicators (1 for event, 0 for censoring) \n", " at each time point (assumed to be sorted).\n", " alpha: significance level (e.g., 0.05 for 95% confidence interval)\n", " method: if not \"log-log\" then the function return the \"plain\" \n", " or \"linear\" confidence intervals\n", "\n", " Returns:\n", " A tuple containing two numpy arrays:\n", " standard error using the given method (log-log by default) and alpha (95% by default).\n", " lower and upper bounds of the confidence interval using the given method and alpha.\n", " \"\"\"\n", "\n", " # Use inverse normal cumulative distribution for z-score\n", " z = norm.ppf(1 - alpha/2) # Two-tailed test\n", "\n", " # first we calculate the term in the sum for the variance;\n", " # di is indeed the array `event_indicators_sorted`\n", " \n", " # we need the number of individual at risk remaining at each time point\n", " n = len(d:=event_indicators_sorted)+1 - np.cumsum(np.ones_like(d))\n", "\n", " # and then compute the cumulative sums of the terms\n", " sum_e = np.cumsum(d / (n * (n-d)))\n", "\n", " # then we compute the cumulative sums of the log of the term for the log-log method\n", " sum_l = np.cumsum(np.log((n - d)/n))\n", " \n", " # if the method is different from \"log-log\", then we return the \"plain\" or \"linear\" SE and CI\n", " if method != \"log-log\":\n", " SE = survival_prob_sorted * np.sqrt(sum_e)\n", " return SE, np.array([survival_prob_sorted - z*SE, survival_prob_sorted + z*SE])\n", " else:\n", " SE = np.sqrt(sum_e / sum_l**2)\n", " return SE, np.array(\n", " [survival_prob_sorted**np.exp(1.96*SE), survival_prob_sorted**np.exp(-1.96*SE)])" ] }, { "cell_type": "code", "execution_count": 64, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[0.33404067 0.33404067 0.21278887 0.11775418 0.04810556 0.04810556\n", " 0.04810556]\n", "[0.97856182 0.97856182 0.91211394 0.81325222 0.68548852 0.68548852\n", " 0.68548852]\n" ] } ], "source": [ "# Example usage (assuming we have our Kaplan-Meier estimates)\n", "survival_prob_sorted, event_indicators_sorted, _ = kaplan_meier(data['YEARS'], data['EVENT'])\n", "\n", "SE, (lower_bound, upper_bound) = greenwood_confidence_interval(\n", " survival_prob_sorted, event_indicators_sorted)\n", "\n", "# lower_bound and upper_bound contain the confidence interval for each time point\n", "print(lower_bound)\n", "print(upper_bound)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Conclusion\n", "\n", "Survival analysis is a branch of statistics that deals with time-to-event data. It's used in a variety of fields, from medical research to machine failure analysis. The `lifelines` library provides a high-level, intuitive API for survival analysis in Python, making it accessible for both statisticians and non-statisticians.\n", "\n", "Here are some key points about `lifelines`:\n", "\n", "- Easy to use: it provides a simple and intuitive interface for fitting survival models, calculating survival probabilities, and visualizing survival data.\n", "- Comprehensive: it includes a variety of survival models, such as the Kaplan-Meier estimator and the Cox proportional hazards model, allowing us to choose the model that best fits our data.\n", "- Handles censoring: One of the key challenges in survival analysis is dealing with censored data. It handles right-censored data natively, making it easier to work with this type of data.\n", "- Confidence intervals: it automatically calculates confidence intervals for survival function estimates, providing a measure of uncertainty around these estimates.\n", "- Plotting capabilities: it includes functions for plotting survival functions and hazard functions, providing a visual way to understand our survival data.\n", "\n", "In this chapter, we've seen how to use lifelines to calculate survival probabilities, create survival tables, and estimate confidence intervals. We've also discussed some of the theory behind survival analysis, including the concept of \"at risk\" individuals and the calculation of survival fractions.\n", "\n", "In conclusion, lifelines is a powerful tool for survival analysis in Python. Whether we're a researcher studying patient survival times, a data scientist predicting customer churn, or an engineer analyzing machine failure times, lifelines has the tools we need to analyze our data and draw meaningful conclusions." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Cheat sheet\n", "\n", "This cheat sheet provides a quick reference for essential code snippets used in this chapter.\n", "\n", "### Survival function\n", "\n", "```python\n", "import lifelines\n", "\n", "# Create a fitter instance\n", "kmf = KaplanMeierFitter()\n", "\n", "# Fit the data to the model\n", "kmf.fit(\n", " durations=durations,\n", " event_observed=event_observed,)\n", "\n", "# Median survival time\n", "kmf.median_survival_time_\n", "\n", "# Survival function\n", "kmf.survival_function_\n", "\n", "# Survival table\n", "from lifelines.utils import survival_table_from_events\n", "\n", "survival_table_from_events(durations, event_observed)\n", "\n", "# Survival events\n", "from lifelines.utils import survival_events_from_table\n", "\n", "# time, events, censored\n", "T, E, W = survival_events_from_table(\n", " data.set_index('YEARS'),\n", " observed_deaths_col='observed',\n", " censored_col='censored')\n", "```\n", "\n", "### Confidence interval\n", "\n", "```python\n", "# Confidence interval\n", "kmf.confidence_interval_\n", "```\n", "\n", "### Kaplan-Meier plot\n", "\n", "```python\n", "# Plotting the survival function\n", "kmf.plot(\n", " show_censors=True,\n", " ci_show=True,\n", " at_risk_counts=True, # add summary tables under the plot\n", ")\n", "```" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Session information\n", "\n", "The output below details all packages and version necessary to reproduce the results in this report." ] }, { "cell_type": "code", "execution_count": 65, "metadata": { "tags": [ "hide-input" ] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Python 3.12.7-------------\n", "\n", "numpy: 1.26.4\n", "pandas: 2.2.2\n", "lifelines: 0.29.0\n", "scipy: 1.14.1\n", "statsmodels: 0.14.2\n", "matplotlib: 3.9.2\n" ] } ], "source": [ "!python --version\n", "print(\"-------------\")\n", "\n", "from importlib.metadata import version\n", "\n", "# List of packages we want to check the version\n", "packages = ['numpy', 'pandas', 'lifelines', 'scipy', 'statsmodels', 'matplotlib']\n", "\n", "# Initialize an empty list to store the versions\n", "versions = []\n", "\n", "# Loop over the packages\n", "for package in packages:\n", " try:\n", " # Get the version of the package\n", " package_version = version(package)\n", " # Append the version to the list\n", " versions.append(package_version)\n", " except Exception: # Use a more general exception for broader compatibility\n", " versions.append('Not installed')\n", "\n", "# Print the versions\n", "for package, version in zip(packages, versions):\n", " print(f'{package}: {version}')" ] } ], "metadata": { "kernelspec": { "display_name": ".env", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.12.7" } }, "nbformat": 4, "nbformat_minor": 2 }